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STABILITY AND TRANSITIONS OF THE SECOND GRADE 

POISEUILLE FLOW 

SAADET OZER AND TAYLAN SENGUL 


Abstract. In this study we consider the stability and transitions for the 
Poiseuille flow of a second grade fluid which is a model for non-Newtonian 
fluids. We restrict our attention to flows in an infinite pipe with circular cross 
section that are independent of the axial coordinate. 

We show that unlike the Newtonian (e = 0) case, in the second grade 
model (e / 0 case), the time independent base flow exhibits transitions as the 
Reynolds number R exceeds the critical threshold R c ~ 4.124e —1 / 4 where e 
is a material constant measuring the relative strength of second order viscous 
effects compared to inertial effects. 

At R = R c , we find that generically the transition is either continuous or 
catastrophic and a small amplitude, time periodic flow with 3-fold azimuthal 
symmetry bifurcates. The time period of the bifurcated solution tends to 
infinity as R tends to R c . Our numerical calculations suggest that for low e 
values, the system prefers a catastrophic transition where the bifurcation is 
subcritical. 

We also find that there is a Reynolds number R^; with R e < Rc such that 
for R < R e , the base flow is globally stable and attracts any initial disturbance 
with at least exponential speed. We show that He ~ 12.87 at e = 0 and He 
approaches R c quickly as e increases. 


1. Introduction 

Certain natural materials manifest some fluid characteristics that can not be 
represented by well-known linear viscous fluid models. Such fluids are generally 
called non-Newtonian fluids. There are several models that have been proposed 
to predict the non-Newtonian behavior of various type of materials. One class of 
fluids which has gained considerable attention in recent years is the fluids of grade 
n mnmmmmm - A great deal of information for these types of fluids can 
be found in 0], Among these fluids, one special subclass associated with second 
order truncations is the so called second-grade fluids. The constitutive equation of 
a second grade fluid is given by the following relation for incompressible fluids: 

t = —pi + pAi + Q 1 A 2 + a 2 Aj, 

where t is the stress tensor, p is the pressure, p is the classical viscosity, or and a 2 
are the material coefficients. Ai and A2 are the first two Rivlin-Ericksen tensors 
defined by 

A, = Vv + Vv r , 

A 2 = A, + AiVv + Vv 1 Ai, 

where v is the velocity field and the overdot represents the material derivative with 
respect to time. This type of constitutive relation was first proposed in [2]. The 
following conditions: 

ol\ + a 2 = 0, p > 0, or > 0, 
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must be satisfied for the second-grade fluid to be entirely consistent with classical 
thermodynamics and the free energy function achieves its minimum in equilibrium 

0 - 

Equation of motion for an incompressible second grade Rivlin Ericksen fluid is 
represented as: 

p(y t + w x v + V^-) = —X/p + pAv + a[Av t + Aw x v + V(v • Av + i|Ai| 2 )], 
V • v = 0. 

where p is the density, a = aq = —< 22 , represents the second order material con¬ 
stant. Subscript t denotes the partial derivative with respect to time, w is the usual 
vorticity vector defined by 

w = V x v. 

We next define the non-dimensional variables: 

v* = 0 p' = P ~U\ x* = i 

Up L L 

where U and L are characteristic velocity and length, respectively. By letting e 
represent the second order non-dimensional material constant which measures the 
relative strength of second order viscous effects compared to inertial effects and 
defining the Reynolds number, 

pU L a 


R = 


e = 


p pL 2 ’ 

the equation of motion, with asterisks omitted, can be expressed as: 


(1) 


X/p = —Av + e(Aw x v + Av t ) — v t — w x v, 
R 


where the characteristic pressure p is defined as: 

P = P+ - e(vAv + A L | 2 ). 

Taking curl of both sides of 0 we can simply write the equation of motion as: 


( 2 ) 


V x [—Av + e(Aw x v + Av t ) — v t — w x v] = 0, 
R 


which is the field equation of incompressible unsteady second grade Rivlin-Ericksen 
fluid independent of the choice of any particular coordinate system. 

Now we restrict our interest to flows in a cylindrical tube and assume that the 
velocity is dependent only on two cross-sectional variables x, y and the time t. The 
incompressiblity of the fluid allows us to introduce a stream function ip such that 
v = (ip y , — ip x , w) where ip = 0(t, x, y) and also w = w(t, x, y ). 

We further take the cross section of the cylinder to be a disk with unit radius 
and consider the no-slip boundary conditions. Then the equations (0 admit the 
following steady state solution 

wo = ^-x 2 -y 2 ), -00 = 0 . 

Here the characteristic velocity has been chosen as U = First considering the 
deviation w' = w — wq and 0' = ip — ipo and then introducing the polar coordinates 
w"(t,r,0) = w'(t, r cos6,rsm0) and 0"(t,r, 9) = ip'(t, r cos 9, r sin 9) and ignoring 
the primes, the equations become 
d 1 

■ 777(1 — eA)w = — A w + Ripe + J(ip, (1 - eA)w), 
at R 

d 1 

— A(eA — 1)0 = ——A 2 ip + eRAw$ + J((l — eA)A0, 0) + eJ(Aw, w), 
at R 


( 3 ) 
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in the interior of the unit disk f 1 where J is the advection operator 

J(f,g) = -Urge - fegr)■ 

r 


The field equations are supplemented with no-slip boundary conditions for the 
velocity field 


( 4 ) 


w = ip 


chp_ 

dr 


0 at r = 1. 


In this paper, our main aim is to investigate the stability and transitions of Q 
subject to ©. We first prove that the system undergoes a dynamic transition at the 
critical Reynolds number R c ss 4.124e -1 / 4 . As R crosses R c the steady flow loses 
its stability, and a transition occurs. If we denote the azimuthal wavenumber of an 
eigenmode by m, then two modes, called critical modes hereafter, with m = 3 
and radial wavenumber 1, become critical at R = R c . Using the language of 
dynamic transition theory m, we can show that the transition is either Type- 
I(continuous) or Type-II(catastrophic). In Type-I transitions, the amplitudes of the 
transition states stay close to the base flow after the transition. Type-II transitions, 
on the other hand, are associated with more complex dynamical behavior, leading 
to metastable states or a local attractor far away from the base flow. 

We show that the type of transition preferred in system © is determined by 
the real part of a complex parameter A which only depends on e. In the generic 
case of nonzero imaginary part of A, there are two possible transition scenarios 
depending on the sign of the real part of A: continuous or catastrophic. In the 
continuous transition scenario, a stable, small amplitude, time periodic flow with 
3-fold azimuthal symmetry bifurcates on R > R c . The time period of the bifurcated 
solution tends to infinity as R approaches R c , a phenomenon known as infinite 
period bifurcation jfb. The dual scenario is the catastrophic transition where the 
bifurcation is subcritical on R < R c and a repeller bifurcates. In the non-generic 
case where the imaginary part of A vanishes, the limit cycle degenerates to a circle 
of steady states. 

The transition number A depends on the system parameter e in a non-trivial way, 
hence it is not possible to find an analytical expression of A as a function of e. So, 
A must be computed numerically for a given e. Physically, the transition number 
can be considered as a measure of net mechanical energy transferred from all modes 
back to the critical modes which in turn modify the base flow. We show that A is 
determined by the nonlinear interactions of the critical modes (m = 3) with all the 
modes having m = 0 and m = 6. Moreover, our numerical computations suggest 
that for low e fluids (e < 1), just a single nonlinear interaction, namely the one 
with m = 0 and radial wavenumber 1 mode, dominates all the rest contributions 
to A. Our numerical experiments with low e, i.e. e < 1, suggest that the real part 
of A is always positive indicating a catastrophic transition at R = R c . 

We also determine the Reynolds number threshold R^ > 0 below which the 
Poiseuille flow is globally stable, attracting all initial conditions with at least expo¬ 
nential convergence in the ffg(ff) norm for the velocity. We find that Re ~ 12.87 
when e = 0. The gap between R^ and R c shrinks to zero quickly as e is increased. 

The paper is organized as follows: In Section 3, the linearized stability of the 
system is studied and the principle of exchange of stabilities is investigated. In 
Section 4, transition theorem of the system is presented with its proof given in 
Section 5. Section 6 is devoted to the energy stability of the system. In Section 7, 
we give a detailed numerical analysis. Finally in Section 8 the conclusions and 
possible extensions of this study are discussed in detail. 
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2. The model and the functional setting 


Throughout Re z, Im z, ~z will denote the real part, imaginary part and the 
complex conjugate of a complex number z. 

Let 

Xi = Hq(Q) x x = L 2 (n) x L 2 (n), 

where fl is the unit disk in R 2 , Hq(Q) and H^fl) denote the usual Sobolev spaces 
and L 2 (fl) is the space of Lebesgue integrable functions. 


For (hi = 


Wi(r,6 ) 
ipi(r,0) 


£ X, i = 1,2, the inner product on A' is defined by 


(5) 


(0i, 02) = 


/*27T 

/ 0 Jo 

12 


(wiW 2 + 0102 )rdrdd, 


M = 

I-e A 

0 

, N = 

’ rA 

H/0 


0 

A(eA - 1) 


eRA<9g 

(N 

<1 


with the norm on X defined as ||0|| 2 = (0,0). We define the linear operators 
M : A'i -> A and N : A0 -lias 

( 6 ) 

and the nonlinear operator H : Xi —> X as 

J (0, (1 — eA)w) „ 

1J(0)= V V , for 0 = € X\. 

J((l — eA)Au>, ip) + eJ(Aw, w) 0 

We will use H to denote both the nonlinear operator as well as the bilinear form 

<7(0/, (1 - eA )wj) 

J((l - eA)Aui/, 0j) + eJ(Awi, wj) 


(7) 


7?(0/,0j) = 


for 0/ = [wi,i/ji] t , cf>j = [wj,i/jj] t . Also we will use the symmetrization of this 
bilinear form 


(8) H s (0/, 0j) = Lf (0/, 0j) + H (0j, 0/). 

Then the equations © and ® can be written in the following abstract form 

(9) M0t = iV0 + Lr(0), 0 e Xi, 

with initial condition 

0(0) = 0o G AT,. 


3. Linear Stability 

To determine the transitions of ©, the first step is to study the eigenvalue 
problem N(f> = /3M0 of the linearized operator. This is equivalent to the problem 

—Aw + R0e = /3(1 — eA)w, 

m “ 

—A 2 0 — eRA-icg = /3(1 — eA)A 0, 

R 

subject to the boundary conditions ®- 

An interesting feature of the eigenvalue problem is the following. 

Lemma 1. Any eigenvalue ft of m with boundary conditions © is real. 
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Proof. Multiplying the first equation of m by Aw and the second equation by ip, 
integrating over the domain Q, we obtain after integration by parts 

(11) (4 + /3e)l|AH| 2 + R [ ipeAwdrdO = -/3||Vw|| 2 , 

R Jn 

and 

(12) - (4 + ^e)||A^|| 2 + eR / A wgipdrdd = fd\\\/ip\\ 2 . 

R Jn 

Let 

■A.\ = e||Aw|| 2 + ||Ai/’|| 2 , 

M 2 = e||Vu;|| 2 + ||VV;|| 2 

A 3 = 2e / A weipdrdd. 

Jn 

Now consider —ex ciD + m which is 

(13) — (— + /3e)Ai + Ri?e(Ma) = PA 2 , 

K 

after integrating by parts. Taking the imaginary part of ED gives 

Im(/3)(eMi + A 2 ) = 0. 

Since {eA\ + A 2 ) > 0, we must have Im(/3) =0. □ 


Now we turn to the problem of determining explicit expressions of the solutions 
of the eigenvalue problem of the linearized operator. Thanks to the periodicity in 
the 9 variable, for to £ Z and j £ Z + , we denote the eigenvectors of ED by 


(14) 


(pm,j{r, 9) — e (pm jfr), — 


Wm,j{r) 


with corresponding eigenvalues fd m j. Let us set the eigenvalues fJ m j for m 0 to 
be ordered so that 

Pm, 1 'll pd m 2 ^ ’ 1 

for each to £ Z. 

Plugging the ansatz ED into ED and omitting j we obtain two ODE’s in the 
r-variable. 

+ fie ) A m Wm + irnRipm = (dw m , 

(15) 


R 

- ^ + Pcj A 2 m ip m + iemRA m w m = -f)A m ip m , 

with boundary conditions 

(16) w m (l) = ip m {l) = = 0 

where 

d 2 1 d to 2 

m dr 2 ^ r dr r 2 

When to = 0, the equations ED become decoupled and we easily find that there 
are two sets of eigenpairs given by 

— ..2 

w o,j( r ) = Maojr), V’ojM = °> A \,j = 


x 0,j 


R(1 + ea qj) 


Woj(r) = 0, ipZjir) = J 0 (a ljr ) - J 0 (a hj ), fd^■ = 


*rj 


R(l + ea? •)’ 
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where otkj is the j th zero of the kth Bessel function </&. In particular, Bo,j < 0 for 
all e and for all R. 

In the / 0 case, as the eigenvalues are real by Lemma |1] the multiplicity 
of each eigenvalue B = Bm,j = B-m,j £ K is generically two with corresponding 
eigenvectors (f>mj and = 4>m,j- 

Solving the first equation of © for -i/> m and plugging it into the second equation 
of (fl5l) yields a sixth order equation 

(17) (A + A m )(p, + A m )A m w m = 0 

where 

n , y/emR - Bm,j —y/emR - Bm,j 

18 J_ -u e /3 — + eB ■ ' 

It is easy to check that A = 0 or /i = 0 yield only trivial solutions to equations © 
and m So we will assume A 7 ^ 0, and /i 7 / 0. 

When to > 0, the general solution of © is 

Wm = c\r m + c 2 Jm(V\r) + c 3 J m (y/jIr) + c 4 r _m + c 5 Y m (\/Xr) + c 6 Y m (^JIr), 


where J m and Y m are the Bessel functions of the first and the second kind respec¬ 
tively. 

The boundedness of the solution and its derivatives at r = 0 implies that c 4 = 
C 5 = cq = 0 and we get the eigensolutions 


(19) 


Wm,j{r) 


V’mj(r-) 


c 4 r m + c 2 J m (VXr) + c 3 J m (y/]ir) 
■> 


if to > 0 
if to < 0 


dir m + d 2 Jm(VXr ) + d 3 J m (y/jdr) 


if TO > 0 
if to < 0 


with di 


■ 

toR 


d 2 = —i^/ec 2 , d 3 = iy/ec 3 . 


The eigenvalues and two of the three coefficients Ci, c 2 , c 3 in © are determined 
by the boundary conditions © which form a linear system for the coefficients Ck- 
This system has a nontrivial solution only when the dispersion relation 



1 dm ( \/A) 

(20) 

B y/emRJ m (y/\) -y/emRJ m (yfjl) 


fim y/emRVXJ^iV A) -y/lmR'/JlJ' rn ( y /Jl) 

is satisfied. 


Using the identity 


ill 

(^) J m (-2’) Jm-\-1 (^) 

Z 

we can show that (OIH) is equivalent to 

(21) 

V^A Jm(V^) t An+l (\/m) 


where J m is the Bessel function of the first kind of order to. 

To compute the critical Reynolds number R c , we set B = 0 in m which, after 
some manipulation, yields 

( 22 ) = 0 . 


Once © is solved for A, the corresponding Reynolds number is obtained by the 
relation A = yj ctoR 2 (from (IT5T) when /3 = 0). We note here that this is the exact 
same equation as the one obtained in ©■ For each in, the equation © has 
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m 

1 

2 

3 

4 

5 

6 

^m, 1 

21.260 

34.877 

51.030 

69.665 

90.739 

114.21 

(A mp/m ) 1 / 2 

4.610 

4.175 

4.124 

4.173 

4.260 

4.36 


Table 1. The smallest positive root A m p 
for m = 1,..., 5. 


of 


and (A 


m.l 


/ m ) 1 / 2 


infinitely many solutions {AmjJ^j where A m ,j increases with j and X m j —> oo as 
j —► oo. Letting R m j = e” 1 / 4 (Awe have /3 m j = 0 when R = R m j. 

We define the critical Reynolds number 


R c = min R m j = min R m .i 

m,j ’ m 


e “ 1/4 min 



so that 8 m,j < 0 if R < R c . There has been recent progress on the properties of 
zeros of K21) . In [I], the estimate 

2 4 (to + 1)(to + 2 )(to + 3 )-\/(m + 4)(m + 5) 2 

\/5 to +15 171,1 

2 4 (m + l)(m + 2)(?n + 3 )(to + 4 )(m + 5) 

< 5m +17 ’ 

on A m y is obtained. 

Using we can show that the upper bound for /m for m = 3 is less 

than its lower bound for all m > 7 which implies that -\/Aiyi/3 < sj A mj i/rn for all 
m > 7. Thus R c is minimized for some m smaller than 7 and hence can be found 
by brute force. Looking at Tabic |TJ we find that the expression above is indeed 
minimized when m = 3 and obtain the relation (fiHl) . i.e. R c = R 3 ,i ~ 4.124e -1 / 4 . 

Defining the left hand side of (12TT) as w(/3, R), the equation (12T1) becomes w(/3, R) = 
0. By the implicit function theorem, this defines 83,1 (R) for R near R c with 
/^ 3 ,i(Rc) = 0 . With the aid of symbolic computation, we can compute 

d&,i , du/dR 0.12^ 

dR ' R=Rc dw/dft l R=R <=>^=^ 3 . 1 0.02 + e > ' 

Thus we have proved the Principal of Exchange of Stabilities which we state below. 


Theorem 1. For e ^ 0, let 


(24) 

Then 


(25) 


Rc = e _1/4 



4.124e" 1/4 . 


{ <0 if R < R c 
= 0 ifR = R c 
>0 if R> R c 
Pm,j(R c ) < 0, if ( m,j) ^ (±3,1). 


Note that Theorem |T] is in contrast to the e = 0 case where the basic flow is 
linearly stable for all Reynolds numbers. This can be seen easily by noting that 
when e = 0 , the inner product of the second equation in m with yields 

(26) ||AV >|| 2 = —/3R|| V-0|| 2 . 

With the Dirichlet boundary conditions, from (BH1) . it follows easily that /3 < 0 if 

V> 7^ 0. 
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For the proof and presentation of Theorem^ we also need to solve the eigenvalue 
problem of the adjoint linear operator which yields adjoint modes orthogonal to the 
eigenmodes of the linear operator. Adjoint problem is obtained by taking the inner 
product of © by <f>* and moving the derivatives via integration by parts onto (f>* by 
making use of the boundary conditions. This yields the following adjoint problem 


N* <\>* = P*M*(j>*, 


where 
M* = M = 


and iv*, cf>* satisfies the same boundary conditions © as w and </>. We denote 

n T 


I-t A 0 


4 A — eRA<9e 


w 


, N* = 

R u 

, </>* = 


0 A(eA - 1) 


Rde -gA 2 


r 


the adjoint eigenvectors by </= 


r. 


and we also have the adjoint 


eigenvalues The reason we introduce the adjoint eigenmodes is to 

make use of the following orthogonality relation 

(27) (<j>m,i,M<t>* ntj ) =0 if (to, i) £ ( n,j). 


4. Dynamic Transitions 


Let us briefly recall here the classification of dynamic transitions and refer to [10] 
for a detailed rigorous discussion. For e ^ 0, as the critical Reynolds number R c 
is crossed, the principle of exchange of stabilities E5l) dictates that the nonlinear 
system always undergoes a dynamic transition, leading to one of the three type of 
transitions, Type-I(continuous), Il(catastrophic) or Ill(random). On R > R c , the 
transition states stay close to the base state for a Type-I transition and leave a local 
neighborhood of the base state for a Type-II transition. For Type-III transitions, a 
local neighborhood of the base state is divided into two open regions with a Type-I 
transition in one region, and a Type-II transition in the other region. Type-II and 
Type-III transitions are associated with more complex dynamics behavior. 

Below we prove that for ©, only two scenarios are possible. In the first sce¬ 
nario, the system exhibits a Type-I (or continuous) transition and a stable attractor 
will bifurcate on R > R c which attracts all sufficiently small disturbances to the 
Poiseuille flow. We prove that this attractor is homeomorphic to the circle S' 1 and 
is generically a periodic orbit. The Figure [Tj shows the stream function of the bi¬ 
furcated time-periodic solution. The dual scenario is that the system exhibits a 
Type-II (or catastrophic) transition. 

The type of transition at R = R c depends on the transition number 


(28) A — 'y ) Aqj + Aqj, 

o =i 

where A m j represent the nonlinear interaction of the critical modes with the mode 
with azimuthal wavenumber to and radial wavenumber j. The formulas for A m j 
are 


A o,j = 


1 


-<F O y(# s (03,l,0Oj),</>3,l) 


= AA- afA* \ ^ 6 b(- g s(</>3,l,<( > 6,i),</ ) 3,l) 

\93,l> M 93,l) 

= 


fio,j {4*0,j ! 

1 


3,1,03,l), <Po ,j)> 


®6,j - — - J2 -3,lj 03 ,i)j 06,j)» 

,ji lvl V6,j> 


( 29 ) 
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Figure 1. The time periodic stream function if per given by (®D 
which rotates in time, clockwise if T > 0 and counterclockwise if 
T < 0. 


To recall the meaning of various terms in (l29l) . /3 n j is the jth eigenvalue of the 
nth azimuthal mode with corresponding eigenvector (j) n j and adjoint eigenvector 
4>n,j- (">') denotes the inner product (©, H denotes the bilinear form Q, H s is its 
symmetrization (©, M is the linear operator defined by ©• 


Theorem 2. If e ^ 0 then the following statements hold true. 

(1) IfRe(A) < 0 then the transition at R = R c is Type-I and an attractor E^ 
bifurcates on R> R c which is homeomorphic to S 1 . Iflm(A) = 0 then Y,r 
is a cycle of steady states. If Im(A) ^ 0 then Y,r is the orbit of a stable 
limit cycle given by 


(30) 


'Ulperifi T, $) 
fjper{ti T, 0) 




w 3 ,i (r) 

1p3,l(r) 


+ o(\/— 


with time period 


(31) 


27rRe(T) 
Im(A)/3 3a ' 


Here w; 3j i and ?/> 3j i are the vertical velocity and stream function of the eigen¬ 
mode of the linear operator with corresponding eigenvalue /3 3 > i- 
(2) If Re (A) > 0 then the transition at R = R c is Type-II and a repeller E^ 
bifurcates on R < R c . Iflm(A) = 0 then E^ is a cycle of steady states. If 
Im(A) 0 then Yir is the orbit of of an unstable limit cycle given by (1301) 
with /3 3j i replaced by —/3 3 j i- 


Remark. In the generic case o/Im(A) ^ 0, Theorem^ guarantees the existence of 
a stable (unstable) bifurcated periodic solution on R > R c (R < R c ). By (1311) . the 
period of the bifurcated solution approaches to infinity as R \. R c (Rf R c ). 


5. Proof Of Theorem [2] 

As is standard in the dynamic transition approach, the proof of Theorem [ 2 ] 
depends on the reduction of the field equations © on to the center manifold. 













10 


OZER AND SENGUL 


Let us denote the (real) eigenfunctions and adjoint eigenfunctions corresponding 
to the critical eigenvalue , 83,1 by 

ei(r , 9) = Re(0 3 j i(r,6»)) e 2 (r,0) = Im(</> 3 i i(r, 9)) 

e i(r,6) = Re(^ 3 > 1 (r, 9)) e 2 (r,9) = Im(c (>3 >:L (r, 9)) 

By the spectral theorem, the spaces X\ and X can be decomposed into the direct 
sum 

X x = E\ ® E 2 , X = E\ ® E 2 , 

where 


E ± = span{ei,e 2 }, 

E 2 = {u&Xi\ {u, e*) = 0 i = 1,2}, 
E 2 = closure of E 2 in X. 


Since M : Xi —> X is an invertible operator, we can define L = M 1 iV and 
G = M _1 7L. Now the abstract equation ([9]) can be written as 

(32) = L(j> + G(4>). 

The linear operator L in (l32l) can be decomposed into 

L = J® C, 

J = L \e 1 ' Ei —> Ei, 

£ = L \e 2 '- E 2 —► E 2 , 

Since the eigenvalues are real, we have Lek = /3 3 j i for k = 1,2. Hence we have 
J = /3 3 j i(R)/ 2 , where J 2 is the 2x2 identity matrix. We know that when J 
is diagonal, we have the following approximation of the center manifold function 
$ : Ei —> E 2 near R ss R c ; see ma- 

(33) - C$(x) = P 2 Gk{x) + o(k), 

The meaning of the terms in the above formula (1331) are as follows. 

a) o(fc) = o(||a;|| fc ) + 0(|/3 3i i(R)|||a;|| fc ) as R R c , x ->• 0, 

b) P 2 : X —> E 2 is the canonical projection, 

c) x is the projection of the solution onto E\, 

(34) x(t, r, 9) = xi(t)ei(r , 9) + x 2 {t)e 2 (r, 9) 

d) Gk denotes the lowest term of the Taylor expansion of G(u) around u = 0. 
In our case G is bilinear and thus k = 2 in and G = Gk ■ 

It is easier to carry out the reduction using complex variables. So we write (1311) 
as 


(35) x{t, r, 9) = z(t)(j> 3i i(r, 9) + z(t)<j> 3A (r, 9) 

where z(t) = ^(xi(t) — ix 2 {t)). Let us expand the center manifold function by 


(36) $= ®n,j(t)<l>n,j(r,9) 

(n,j)#(± 3,1) 


Plugging the above expansion into the center manifold approximation formula (1331) . 
taking inner product with Mq i>* ■ and using the orthogonality (1271) we have 


(37) 


^n.j — 


1 


/3n,j{4>n,j, 


{H(x), <j>n,i) + o(2). 


Since H is bilinear, 

(38) _ __ 

H(x) = H(z4 > 3 ,1 + z(/) 3t i) = z 2 H{cf> 3 } i, <f> 3 ,\) + zzH s {(j) 3 ^i,^ 3j i) + z 2 H(4> 3} i, </> 3 j i), 
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with the operator H s defined by Thanks to the orthogonality e m6 e lrn6 dd = 
2 n'S nm , we have 

(39) = °i if TOi + TO 2 ^ m 3 

With <f> 3j i = this implies 

(40) (H(x),r„, j ) = 0 if ni { 0 ,- 6 ,6}. 

According to (TSU) . (1M1) and (IdTl) . 

OO 

( 41 ) <&(t) = ^2 ®0,j(t)<f>0,j + + 0 ( 2 ) 

7=1 

That is the center manifold is o(2) in eigendirections whose azimuthal wavenumber 
is not 0, 6 or —6. The equation (1551) implies that 

(H{x),(j> oj) = zz{H s (tj> 3tl ,</> 3il ),0Zj), 

(42) (H(x), r 6j ) = z 2 (H(^ i, 0 3i i), r 6j ), 

By (1421) and (1571) . we get the coefficients of the center manifold in m 

= ZZ® o,j + o( 2), 

(43) =z 2 $ 6j +o(2), 

^—6 j = 4>6, j, 

where and ‘4 , 6,j are given by (EH). 

As the dynamics of the system is enslaved to the center manifold for small 
initial data and for Reynolds numbers close to the critical Reynolds number R c , 
it is sufficient to investigate the dynamics of the main equation E on the center 
manifold. For this reason we take 

<t>{t) = x(t) + 4>(f), 

in © to obtain 

(2% dtz _ _ 

(44) — M(j) 3 i + — Mc/) 3 i = zNfai + zNfo ! + H(x + $). 

at at 

To project the above equation onto the center-unstable space we take inner 
product of (Hill with <j)^ 1 and use 

(Mfo2,ct>l 1 )= 0 , 

and 

N(f) 3 , i = P 3 ,iM(j) 3t i, Ncf> 3} i = 

to get the following reduced equation of ©■ 
dz 1 

(45) -=/3 3 A (R)z + - - ——{H(x + *),p 3 tl ). 

at (</> 3 ,i, M(p 31 } 

The reduced equation (l45ll describes the transitions of the full nonlinear system 
for R near R c and small initial data. At this stage, the nonlinear term in (li5l) is 
too complicated to explicitly describe the transition. Thus we need to determine 
the lowest order expansion in z of the nonlinear term ( H(x + 4>),0g 2 ). By the 
bilinearity of H , 

(46) (H(x + $), <^,i> = (H(x),r 3 ,i) + (H s (x, $), r 3 ,i) + 4>h). 
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The first term in (1551) vanish by (1551) and the last term in (1461) is o(3) as ff(<I>) = o(3) 
since $ = 0(2) and H is bilinear. Thus (1551) becomes 

(47) (H(x + <&), = (H s (x, $), 0^> + o(3). 

Using the expression (1551) for x, we can rewrite (1571) as 

(48) (H(x + $), 0^) = z(H B (</> 3 ,i, $), 0^> + z(H s (^, $), 0^) + o(3). 

Now we use the expansion (1411) of $ in (1551) and the orthogonality relations 


(^ s fe,i^n,j)^3,i) = 0 if n ^0, 

and 

(H s (<j> 3 ,i,<l>n,j),<t>t,i) = 0 if n 7^6, 
which follow from (1551) to arrive at 
(49) 

OO 

(Zf(x + $), 03 !) = ^2 Z^0,j(Hs(h,l,ct>0,j), 03,l) +^6j(-ff S (03,l, 06j), 03,l) + o(3). 
2=1 


Defining the coefficient A by (1551) and making use of (1551) and (1551) , we write down 
the approximate equation of © as 

rl 7 

(50) ^=/3 3>1 (R)z + A\z\ 2 z + o(3). 

To finalize the proof, there remains to analyze the stability of the zero solution 
of (1551) for small initial data. In polar coordinates z(t) = \z\e 1J , (1551) is equivalent 
to 

(5!) ^ = /? 3 ,i(R)M +Re(A)| 2;| 3 + o(|z| 3 ), ^ =lm(A)\z\ 2 + o(\z\ 3 ). 

For R > R c as ^ 3 p > 0, it is clear from m that z = 0 is unstable if Re(A) > 0 
and is locally stable if Re(7l) < 0. In the latter case, the bifurcated solution is 


z(t) = 


' -fo,i(R) 

Re(A) 


exp —i 


. Im(A) 
1 Re(yl) 


&,l(R)t 


Thus to determine the stability of the bifurcated state as R crosses the critical 
Reynolds number R c , we need to compute the sign of the real part of A. 

The details of the assertions in the proof of Theorem [5] follow from the attractor 
bifurcation theorem in m ■ That finishes the proof. 


6. Energy Stability 

In this section we study the energy stability of the equations m which is related 
to at least exponential decay of solutions to the base flow. We refer to [15] for a 
multitude of applications of this theory. 

For /, g , h in 77g(U), the following two properties of J follows from integrating 
by parts 

( J(f , g), h) = 0, if /, g and h are linearly dependent. 

Taking inner product of the first equation in J5J with w and the second equation 
in ([5]) with 0 and using the property (1551) . we have 

(53) ^(l|w|| 2 + e||Vu;|| 2 ) = ~^|| Vw;|| 2 + R(0e, w) + e{J(Aw, 0), w), 

(54) i|(||V0|| 2 + e ||A0|| 2 ) = -^||A0|| 2 + eR(Aw e ,0) +e(J(Aw,w),0). 
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Adding equations (l5ifl) and (15H) and using (15^1) once again, we arrive at 
(55) 


“£(*) = - ^Mt) + Kl2(t), 


where 


\w\\~ + e|| Vw|| 


£ = IMI 2 - 

II = 11V ZD 11 2 
Z 2 = (4>e, w — eAw). 


2 + IIVV’II 2 

||AV.|| 2 


-l|AV’|| 5 


Letting 

(56) 

we have by (l55l) 

(57) 


1 2a 

— 5 - = max —, 

R E Xx\{0} 2i 


4 £ < -2R(4 - -4)1 1 . 

dt ~ V R 2 R 2 ' 


Since Xi > 0 and I 2 = 0 whenever ipg = 0, R^ must be nonnegative. 

Since w £ Hq(£1) and V0 £ ifg(fi), by the Poincare inequality, |Vw ;| 2 > 771 |zc | 2 
and |Ai /.| 2 > 771 1Vz/^>| 2 , where 171 ss 5.78 is the first eigenvalue of negative Laplacian 
on f2. Thus we have 

Vi c 


(58) 

Now let 


2i > 


1 + em 


cr = 


2R?7i , 1 
1 + erji R 2 


R1 


-), 


and suppose that R < R^;. Then cr > 0 and by <53 and (EH|), 

< —c R £(t). 

Hence the Gronwall’s inequality implies 

£{t ) < e _CRt £( 0 ). 

In particular, for R < R^, cr > 0 and any initial disturbance in X\ will decay to 
zero implying the unconditional stability of the basic steady state solution. 

Using the variational methods to maximize the quantity in (1561) . we find the 
resulting Euler-Lagrange equations as 

>2 


(59) 


R 

A w + —(1 - eA )if>g = 0 , 

r>2 

A 2 ijj + — (1 — eA )wg = 0. 


Considering (1591) as an eigenvalue problem with R playing the role of the eigenvalue, 
R_e is just the smallest positive eigenvalue. To solve <59|), we plug the ansatz 


AmO . 


Wm(r) and ^ = e irne 'ijj rn (r) into (EH which yields 


A„ 


(60) 


.rriR 


. toR' 


(1 - eA m )V> m = 0 


A m^Pm T 7 (1 eA m )w m 0 , 


W 1 "h 7 dr — 73 - • Taking A m of the second equation above and using 


where A m = d 1 1 d 
the first equation, we obtain 

(61) 


p(A m )VVn = 0 
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where 

»(«) = «“ + 

Let £ 1 , £2 and £3 be the three roots of p. As the discriminant of p is negative, one 
root is real and the others are complex conjugate. The factorization of the operator 
in (ELD gives 

(62) (A m - &)(A m - 6 )(A m - £ 3 )^m = 0. 

The general solution of (f(x2j) is 

3 

'Ipm = ^ ^ CkIm ( \J ) T &kKm (\/) 5 

k =1 

where Im and K m are the modified Bessel functions. The boundedness of the 
solution at r = 0 necessitates c.k =0 for k = 1, 2,3. Thus 

3 

'tpm — ^ ^ Cklmjy/fc.k'r )• 

k =1 

Now applying the operator 1 — eA m to the second equation in (fGOl) and using the 
first equation of (EDD. we obtain p(A m )w m = 0, i.e. the same equation ED, this 
time for w m . Hence 

3 

w m — 'y ) d/;; I rn (\J r ). 

fc =1 

The first equation in (1501) gives the relation d/ c = — 1 ) c fe between Ck s 

and dk s. Now the boundary conditions w m (l) = i/jm{ 1) = V’m(l) = 0 constitute a 
homogeneous system of three linear equations for the coefficients Cfc’s. The existence 
of nontrivial solutions is then equivalent to the vanishing determinant of this system 
which after some manipulation becomes 


(63) 


Im(V£i) Im(V&) 
VTMVTi) V&im(V&) 


Im(VTs) 

VTsiUVTs) 


= o. 


For fixed m and e, the equation (1631) has infinitely many solutions R 
jeZ.|_. Letting 


Hmj (c), 


(64) 


Rm — min Rm,? - , 
J€ Z+ 


the critical Reynolds number is given by 


(65) 


Rb = min R m . 


We present the numerical computations of Rb in the next section. 


7. Numerical Results and Discussion 

7.1. Transitions at the critical Reynolds number R c . As shown by Theo¬ 
rem [2l the transition at the critical Reynolds number R c ss 4.124e -1 / 4 is dictated 
by the real part of A given by (051) . 

Let us define 

A= lim A n , 

N—> OO 

where A N is the series in (051) truncated at N, i.e. A N = ,j + A 6j and 

A m ,j represents the nonlinear interaction of the critical modes and the mode with 
azimuthal wavenumber m and radial wavenumber j given by m- 
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Figure 2. The plots of Re(A JV ) for different e values. All A N, s 
are scaled so that |Re(A 1 )| = 1. 


A symbolic computation software is used to compute A N . We present our numer¬ 
ical computations of A N in Figure [2] for e = 1, ICR 1 , ICR 2 , ICR 3 and 1 < N < 10. 
The imaginary part of A is nonzero and we are only interested in the sign of the real 
part of A to determine the type of transition according to Theorem [2] To simplify 
the presentation, we scale all A^’s so that |Re(A 1 )| = 1. The plots in Figure [2] 
suggest that the convergence of the truncations A N —> A is rapid for small e but 
a higher order truncation (larger N) is necessary to accurately resolve A for larger 
e. For e < 10 _1 , even A 1 is a good approximation to determine the sign of A. For 
example, the relative error for approximating A with A 1 is approximately %2 for 
e = ICR 3 and increases to approximately %18 for e = 1. 

We also measure the relative strength of the nonlinear interactions, i.e. the ratio 


( 66 ) 


B N _ £ j= l Rc (^6j) 

EJLiI MAoj)" 


in Figure [3] 

It is seen from Figure [3] that the contribution from the modes with m = 0 
dominates when e is low. But as e increases, the contribution from modes with 
m = 6 start to become significant. For example, for e = ICR 3 , B N approaches 
8 x ICR 5 , for e = ICR 2 , B N approaches —2 x ICR 3 and for e = ICR 1 , B N approaches 
—8 x ICR 3 . In particular, for low e, we have A « A 1 ss Aop. 

More significantly, our numerical results presented in Figure [2] show that the 
real part of A is positive for e = ICR 3 , ICR 2 , ICR 1 ,1, meaning that the transition 
is catastrophic by Theorem [2] Thus the system moves to a flow regime away from 
the base Poiseuille flow and the system exhibits complex dynamical behavior for 
R > R c . 


7.2. Determination of Energy Stability Threshold Rg. With a standard 
numerics package, R m (e) in (1641) can be computed for given m and e. Then by 
(l65l) . R^ is computed by taking minimum in (1651) over all (computed) R m . In 
Table [2] it is shown that for e = ICR 4 , and e = ICR 3 , R# is obtained for m = 1 
while for e = ICR 2 and e = 2 x ICR 2 , R_e is obtained for m = 2. In Figure[5j we plot 
R-m (m = 1,2,3) for 0 < e< 5 x ICR 2 . We see that the curves Ri and R 2 intersect 
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Figure 3. B N , defined by (1M>1) . measures the relative strength of 
the nonlinear interactions of the critical modes with m = 6 modes 
to the m = 0 modes. 

approximately at e = 0.009 while R 2 and R 3 intersect approximately at e = 0.024. 
As e is increased, the value of m for which R^ is minimized also increases. For 
higher values of to, the roots of the determinant in (16311 becomes increasingly hard 
to find. 

In Table [2] the last column gives the value of R c , the linear instability threshold, 
computed by (El) . Note that the interval [Rb,R c ] consists of Reynolds numbers 
for which the base flow is either not globally stable or globally stable but not not 
exponentially attracting. We plot the Re and R c data from Table [2] in Figure [4] 
which shows that this interval shrinks rapidly as e is increased. 


6 

Rr 

r 2 

r 3 

R4 

r 5 

Re 

R c 

0 

12.87 

13.49 

14.84 

16.37 

17.95 

12.87 

OO 

10" 4 

12.86 

13.47 

14.81 

16.32 

17.88 

12.86 

42.4 

io~ 3 

12.77 

13.31 

14.54 

15.90 

17.29 

12.77 

23.84 

10 -2 

11.99 

11.95 

12.44 

12.95 

13.39 

11.95 

13.40 

2 x 10“ 2 

11.26 

10.83 

10.91 

11.03 

11.13 

10.83 

11.27 


Table 2. R m denotes the first positive root of 


and Re is the 


minimum of R„ 
threshold. 


taken over all to. R c is the linear instability 


8. Concluding Remarks 

In this work, we considered both the energy stability and transitions of the 
Poiseuille flow of a second grade fluid in an infinite circular pipe with the restriction 
that the flow is independent of the axial variable z. 

We show that unlike the Newtonian (e = 0) case, in the second grade model 
(e / 0 case), the time independent base flow exhibits transitions as the Reynolds 
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e -10 -2 

Figure 4. and R c curves in the e — R plane. 



Figure 5. The plot of R m vs e for to = 1,2,3. 


number R exceeds the critical threshold R c ~ 4.124e -1 / 4 where e is a material con¬ 
stant measuring the relative strength of second order viscous effects compared to 
inertial effects. At R = R c we prove that a transition takes place and that the type 
of the transition depends on a complex number A. In particular depending on A 1 
generically, either a continuous transition to a periodic solution or a catastrophic 
transition occurs where the bifurcation is subcritical. The time period of the peri¬ 
odic solution approaches to infinity as R approaches R c , a phenomenon known as 
infinite period bifurcation. 

We show that the number A = ]Cy=i -^o ,j + ,j where A m j denotes the nonlin¬ 
ear interaction of the two critical modes with the mode having azimuthal wavenum¬ 
ber to and radial wavenumber j. Our numerical results suggest that for low e 
(e < 1), A is approximated well by Aq : i. That is, a single nonlinear interaction 
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between the critical modes and the mode with azimuthal wavenumber 0 and ra¬ 
dial wavenumber 1 dominates all the rest of interactions and hence determines the 
transition. 

Our numerical results suggest that a catastrophic transition is preferred for low 
e values (e < 1). This means, an unstable time periodic solution bifurcates on 
R < R c . On R > R c , the system has either metastable states or a local attractor 
far away from the base flow and a more complex dynamics emerges. 

We also show that for R < R^ with Rg < R c , the Poiseuille flow is at least 
exponentially globally stable in the Hq(CI) norm for the velocity. We find that 
R^ w 12.87 when e = 0 and the gap between R^ and R c diminishes quickly as e is 
increased. 

There are several directions in which this work can further be extended. 

First , in this work we consider a pipe with a circular cross section of unit radius 
and find that the first two critical modes have azimuthal wavenumber equal to 3. 
Increasing (or decreasing) the radius of the cross section will also increase (decrease) 
the azimuthal wavenumber of the critical modes. The analysis in this case would 
be similar to our presentation except in the case where four critical modes with 
azimuthal wave numbers m and m + 1 become critical. In the case of four critical 
modes, more complex patterns will emerge due to the cross interaction of the critical 
modes [Eli. An analysis in the light of |14| is required to determine transitions 
in this case of higher multiplicity criticality. 

Second, for the Reynolds number region between R# and R c , there may be 
regions where the base flow is either globally stable but not exponentially attractive 
or regions where the domain of attraction of the base flow is not the whole space. A 
conditional energy stability analysis is required to resolve these Reynolds number 
regimes m- 

Third , the results we proved in this work for the second grade fluids can also be 
extended to fluids of higher grades and to other types of shear flows. 

Fourth , in this work we restricted attention to 2D flows. In the expense of 
complicating computations and results, a similar analysis could be considered for 
3D flows which depend also on the axial variable z. 
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